{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" }, "colab": { "name": "302_4th Order Runge Kutta.ipynb", "provenance": [], "include_colab_link": true } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "2CkrgqciiwHJ" }, "source": [ "# Example 4th order Runge Kutta\n", "\n", "The general form of the population growth differential equation\n", "\\begin{equation} y^{'}=t-y, \\ \\ (0 \\leq t \\leq 2) \\end{equation}\n", "with the initial condition\n", "\\begin{equation}y(0)=1,\\end{equation}\n", "Has the exact soulation. \\begin{equation} y= 2e^{-t}+t-1.\\end{equation}\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "Bn8ulRQfiwHK" }, "source": [ "#### Setting up Libraries" ] }, { "cell_type": "code", "metadata": { "id": "13OrAvU_iwHL" }, "source": [ "import numpy as np\n", "import math \n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt # side-stepping mpl backend\n", "import matplotlib.gridspec as gridspec # subplots\n", "import warnings\n", "\n", "warnings.filterwarnings(\"ignore\")\n" ], "execution_count": 1, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "2WMllWf6iwHR" }, "source": [ "## Defining the function\n", "\\begin{equation}f(t,y)=t-y\\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "--2I3VnQiwHS" }, "source": [ "def myfun_ty(t,y):\n", " return t-y" ], "execution_count": 2, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "O-O26YVeiwHV" }, "source": [ "## Initial Setup\n", "Defining the step size $h$ from the interval range $a\\leq t \\leq b$ and number of steps $N$\n", "\\begin{equation}h=\\frac{b-a}{h}.\\end{equation}\n", "This gives the discrete time steps,\n", "\\begin{equation}t_{i}=t_0+ih,\\end{equation}\n", "where $t_0=a$." ] }, { "cell_type": "code", "metadata": { "id": "Cu78WjZziwHW" }, "source": [ "# Start and end of interval\n", "b=2\n", "a=0\n", "# Step size\n", "N=4\n", "h=(b-a)/(N)\n", "t=np.arange(a,b+h,h)" ], "execution_count": 3, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "4aeaLmHSiwHY" }, "source": [ "## Setting up the initial conditions of the equation\n", "\\begin{equation}w_0=IC\\end{equation}\n" ] }, { "cell_type": "code", "metadata": { "id": "p4nV_AdhiwHZ" }, "source": [ "# Initial Condition\n", "IC=1\n", "w=np.zeros(N+1)\n", "y=(IC+1)*np.exp(-t)+t-1#np.zeros(N+1)\n", "w[0]=IC" ], "execution_count": 4, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "TTLxIu5BiwHb" }, "source": [ "## 4th Order Runge Kutta \n", "\\begin{equation}k_1=f(t,y),\\end{equation}\n", "\\begin{equation}k_2=f(t+\\frac{h}{2},y+\\frac{h}{2}k_2),\\end{equation}\n", "\\begin{equation}k_3=f(t+\\frac{h}{2},y+\\frac{h}{2}k_2),\\end{equation}\n", "\\begin{equation}k_4=f(t+\\frac{h}{2},y+\\frac{h}{2}k_3),\\end{equation}\n", "\\begin{equation}w_{i+1}=w_{i}+\\frac{h}{6}(k_1+2k_2+2k_3+k_4).\\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "BVbEmJiiiwHb" }, "source": [ "for k in range (0,N):\n", " k1=myfun_ty(t[k],w[k])\n", " k2=myfun_ty(t[k]+h/2,w[k]+h/2*k1)\n", " k3=myfun_ty(t[k]+h/2,w[k]+h/2*k2)\n", " k4=myfun_ty(t[k]+h,w[k]+h*k3)\n", " w[k+1]=w[k]+h/6*(k1+2*k2+2*k3+k4)" ], "execution_count": 5, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "HoEhZkVMiwHd" }, "source": [ "## Plotting Results" ] }, { "cell_type": "code", "metadata": { "id": "NqpMeCMOiwHd", "colab": { "base_uri": "https://localhost:8080/", "height": 302 }, "outputId": "32f20564-470c-46e3-9315-45fd44b9b533" }, "source": [ "fig = plt.figure(figsize=(10,4))\n", "# --- left hand plot\n", "ax = fig.add_subplot(1,3,1)\n", "plt.plot(t,w, '--',color='green')\n", "#ax.legend(loc='best')\n", "plt.title('Numerical Solution h=%s'%(h))\n", "\n", "ax = fig.add_subplot(1,3,2)\n", "plt.plot(t,y,color='black')\n", "plt.title('Exact Solution ')\n", "\n", "ax = fig.add_subplot(1,3,3)\n", "plt.plot(t,y-w, 'o',color='red')\n", "plt.title('Error')\n", "# --- title, explanatory text and save\n", "fig.suptitle(r\"$y'=t-y, y(0)=%s$\"%(IC), fontsize=20)\n", "plt.tight_layout()\n", "plt.subplots_adjust(top=0.85) " ], "execution_count": 6, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAEdCAYAAAARsJF3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd5xU1f3/8ddnl95BEBBYsQBiQTQrGHtBRMXekK4CgTUxCRh7C4qRRE2+/gwgKtEdUFFERUURKzYioLhRUSnSe+8ssOf3x72bjMtsg9k9OzPv5+Oxj5m598y571n2Mp+5c+655pxDREREREQCab4DiIiIiIhUJCqQRURERESiqEAWEREREYmiAllEREREJIoKZBERERGRKCqQRURERESiqEAWEREREYmiAllEypWZfW5mf/adQ0REpDAqkEWk3JhZa+DXwFLfWVKZmT1nZqvNrOYB9vMrM3Nm1i9e2UREKgIVyCJSns4H1gNjy3vDZvbHsJjrXt7brkjM7CSgF/Cwc25bjPXNzWyMmS03s11mttDM/mFm9Qu2dc7NAl4DHjCzWmWffl9mdpWZ/T8z+8TMNof/xuX+9yUiyUUFsoiUp/OBp5xzOzxs+1fh7SwP265IhgGbgZEFV5jZEQS/n+uBL4G/AwuA3wNfmNlBMfr7C9AEuLmsAhfjbuC3QHtgmacMIpJkVCCLSLkwsyrAacA/PUX4FbAV+MnT9r0Lh7h0Al4q5EPKCOBg4Gbn3GXOududc+cQFMptCIrrX3DOfQn8APzGzHy8p/wRaA3UAQZ52L6IJCEVyCKyX8xscPh19pBC1rcJv6KfFi46DZjqnFtSfinBzB42MwccBdQC8sLczsx6lcP2S/t7Kst+bwAMGB+j/RFAZ2Ah+36IuQ/YBvQqZNzyi0AGcF5pXkM8OOc+dM7Ndc658t62iCSvSr4DiEjC+iy8PbmQ9f8PSCf4+hugA/B/ZR0qhq+A54A+wOfA1Kh1H5XD9kv7eyrLfjsBe4HpMdqfHd6+65zLi17hnNtiZp8RFNAnA+8XkuU8YEqJ0ouIVGAqkEVkf30F7AA6FlxhZlcTFEuPO+dyAJxzD5ekUzP7A1CvFDlmO+deK2ylc+4lM6tHUCA/55wbXYq+46FUv6ey6jc88tsemBPr5DyCIRRQ+BCUuQQFcmv2LZBnhLdnFBU43v+2IiJlRQWyiOwX59xuM5sBnGFmTZ1zK+C/hdhjwGrg3v3o+g/AoaVo/xzBTApFOTG8/Wo/8hyQsvo97Ue/zQiOKK8opMu64e2mQtbnL9+nwHXObTKznQTDLIpSFv+2IiJxpzHIInIg8r9a/3XUsnuB5sBtzrnCiq1COedaOuesFD99S9DticBu4D/FNQynNXOl+CnJlGJx/z3tR7/5M1Bs2M9tFWc90LCoBmX0bysiEnc6giwiByK/QOsITDSzowhmFfiC4Oifd2ZWCTgO+N45t6sET5kP7CzFJpaXoE1Z/Z5K02/+rBXVCukrv5iuW8j6/OUbC1lfPWobIiIJTQWyiByIzwHH/04Ue4Lga/yb9ndWgTIYp3o0QVFYouEVzrlzS7Htkor772k/+l0d3saayxjgx/C2dSHrW4W3+4xRDqd3qwf8XFRYjUEWkUShAllE9ptzboOZzQF+FV6h7lxgpHPu6wPoNt7jVNuHtweS6YCU0e+ptP2uANbwv5PxCvowvO1sZmnRM1mYWW3gVGA7sWfAaEMwfdzsYiJrDLKIJASNQRaRA/UpUBN4ElgL3HUgnZXBONX8I6abDyRXHJT492Rmz4bjm/vGq9/wiPI0oKGZHRlj/XzgXaAlcFOB1X8OtxEpZAaM/CPYH8ZYF70NjUEWkYSgI8gicqA+AwYQXITjj865sjoJbH/lX1p6mJkdS3DBi++ccy+Xc47S/J7yD17siXO/rwBXElzye16M9VkEwzYeN7NzgTkE45vPJhhaUVhR35lgfuXXS5A3rszsMuCy8GGT8PbXZvZseH+tc+6W8s4lIonNdPEhETkQZnY6wZHJGUDHinhFMzP7LfA7gq/3qwIPOecO6Ej3fmQo8e/JzL4GjgAOLe4DRyn7rQIsARY65/aZPzls0wIYCnQhOPq+AngV+HOsLGZWF1gJTHHOXVZwfVkzs/sJrvRXmEXOuZblk0ZEkoUKZBE5IGY2CbgIONk5N6O49qmqpL+n8KIm64BHnXO3xqvfqPZ3AA8BJx7oGOiwv98BjwOnO+c+PdD+REQqAo1BFpH9Fp4YdjHBiWEqjgtRyt/T6QRzNj8W537z/R1YTHCU+ICYWXXgDuAVFccikkx0BFlESsXMMoDuBEMAehNcgriDc26712AVTFn9nuLRr5mdQTCu+JFCTroraT9tgWuBZ51zC/e3HxGRikYFsoiUipkNIJgxYSMwFfiDc64kF8tIKWX1e9LvX0Sk7KlAFhERERGJojHIZcTMTjezH4tvWWw/C82sUzwyxatfM8sws61mll4GuVysOVpFpPQq6n4uIlLRJWyBHP7Hv9rMakYt62dmH3mM9V/OuU+cc4VdsSouzKy5mb1iZmvNbJOZfVvCCwuUdju/eJN1zi12ztVyzu2N97bixcy6m9kiM9tmZq+ZWYMi2rqw3dbw5+nyzCoVR/i3viPqb2GrmT1Rhts7y8yWFtNG+7lIAirv/08kvhL9QiHpwO8JpiyqMMysknOuJBP8H6gI8A3B3K67gOP430T5KcvMjiEYo3kR8BUwGhgBdCviacc752JdOEFSz8XOufd8h4ii/VwkcRX7/0msmsHM0kvz4bS07aV4CXsEOfQ34JZw3tBfMLOW4ZHBSlHLPjKzfuH9vmb2mZn93cw2mtkCMzslXL4kPDrdJ+q5Vc3sETNbbGarzGxUOMXRf48CmdltZrYS+FfBI0Nm1sLMJprZGjNbl/8p0syOMLMPwmVrzWxcrNdTiJMIzh7f5pzb45z72jn3dtQ2LzGz78LX91F4xvk+LLis7YNRj/+b3cwiQAbwRvjp99aCv1szO8TMJpnZejObZ2b9o/q638xeMrNsM9sS5sks5nV1MrO5Ye5/mpmV8PeRrwfwhnNumnNuK3APcIWZ1S5lPyL/ZWYjzeyVqMfDzex9C9Q3szfD/XtDeL95VNsGZvYvM1sern/Ngm+/3gYOiTq6dEiMTSfrfi6SkgrUH+uA+8P9c6SZTTazbcDZZtY23Kc3hvvUJVF97NPe2wtKUoleIM8EPgL29zKiHYEcgqtFPQ+8SPBmdCTQE3jCzGqFbR8GWgPtw/XNgHuj+moCNCA4yjMgeiMWjOF7E1gEtAyf+2L+auAvwCFAW6AFcH8J808H/mlm3SyY+il6m62BF4A/AI2AyQRvflVK2DcAzrleBHOmXhx+3frXGM1eBJaGr+Eq4CEzOydq/SVhm3rAJKC4r5i6Evw7tAOuIbgsLmZ2WvgfRWE/p4XPP4bgiFv+a5gP5BL8+xVmmpmtDD/EtCwmn6SmIcBx4Zvb6cCNQJ/wynVpwL8I9v8MYAe//DuPADUI/jYPBv4eTq92AbA83LdqFTIbRbLu5yKprCOwAGgMDAuXdQ/v1wb+DbwBvEvwf8bvgHFmFj10M7q95iGPs0QvkCEoUn9nZo3247k/O+f+FX4tMZ6gOB3qnNvlnHuXoKg6MjyCOQD4o3NuvXNuC8Gwjuiv7POA+8Ln7iiwnQ4Ebyp/Co8C7cyfVN85N885NzV83hqCiwOcWcL8VwOfEBwh/dnMZpvZSeG6a4G3wr53A48A1YFTSvH7KZYFl6U9FbgtfF2zgacJ5mfN96lzbnL4e44AxxfT7cPOuY3OucXAhwQfSnDOfeqcq1fET/5/ELWATQX63ETwn0gsZxJ8cDkKWA68aVHfPEjKea3AB6/+AOE8w70I9tGxwO+cc0vDdeucc68457aH/z8MI9yPzawpQSE80Dm3wTm32zn3cSnyJOt+LpIKYv5/QvDB+P+F3wrl1wyvO+c+c87lEbzv1SJ4P8x1zn1AcKDtuqi+/9veObez/F5Sakj4Atk59y3BH83t+/H0VVH3d4T9FVxWi+DITA1gVv4fOfBOuDzfmiL+QFsAi2KNSzazxmb2opktM7PNBG+8DUsSPnyzvd05dwzBp9DZBDujERTki6La5gFLCI5ex9MhQP6HhnyLCmxnZdT97UC1YgrQgu1rFdawEFuBOgWW1QG2xGhLOBQj1zm3kWBM+2EER/MlNV1W4IPXU/krnHP/JjjqY8BL+cvNrIaZPWnBiaGbgWlAvfDboxYE+8iG/QmTxPu5SCoo7P+TJTHaRi87BFgS7tP5Cu5zsfqQOEn4Ajl0H9CfX/7h5F8dqkbUsv09sWUtQbF8TNQfeV3nXHThVtSE0kuAjELeLB4Kn3ucc64OwdCO0o65xTm3luDo0SEEQz2WE3zdC0D4ZtoCWBbj6dso+vdU1GtbDjSwX47vzShkOwfEgqnzthbxc3rY9Duijl6Z2eFAVeCnEm7KsR//BpL8zOwmgr+l5cCtUauGAG2AjuF+fEb+Uwj2/wYW+9yCUk1Enwr7uUiKiLW/RS9bDrQws+g6reA+pwtZlKGkKJDD2QfGAzdHLVtD8IfU08zSzewGgkuz7k//ecBTwN/N7GAAM2tmZueXsIsvgRXAw2ZW08yqmdmp4braBEc8N5lZM+BPJc1lwUlCx5pZpfCNaxAwzzm3juDo1kVmdq6ZVSZ4A98FfB6jq9nAhRacSNSEYDxjtFXA4bEyOOeWhH3+JXxd7QjGZo4t6esoKRdMnVeriJ9PwqbjgIvDgromMBSYWODoFxDMeGFm7cO/kVrAowR/N3PinV8SWzje90GCD7G9gFvNrH24ujbBh+iNFkwpeF/+85xzKwhOxhthwcl8lS241DME+9ZBZla3iO2m1H4uIkAwBnk7wf8zlc3sLOBi/nf+kpSxpCiQQ0OBmgWW9ScoONcRnBwT602jpG4D5gHTw69Q3yM4YlSscEzexQQn9y0mONHl2nD1n4ETCcbIvgVMLEWmGsCrBJecXUBwJOmScJs/EryR/z+CI+AXE5yAkxujn/xppBYSnBAwvsD6vwB3h8NLYp0QeR3BGN7lYZ77nMdpspxz3wEDCQrl1QTFS1b+ejN728zuDB82Jni9mwl+hy2BruF4TklN+TM55P+8Gn77MxYY7pz7xjk3F7gTiJhZVeAfBGN/1xKcVPdOgT57AbuBHwj+Jv8A4Jz7geAkuwXh/hVrFgvt5yKJa5//T0rypHAfvpjg/IW1BFOV9g7/z5ByoEtNi4iIiIhESaYjyCIiIiIiB0wFsoiIiIhIFBXIIiIiIiJRVCCLiIiIiETxNol7w4YNXcuWLX1tXqRMzZo1a61zbn+u7pjQtF9LMkvF/Vr7tCSzovZpbwVyy5YtmTlzpq/Ni5QpM1tUfKvko/1aklkq7tfapyWZFbVPa4iFiIiIiEgUFcgiIiIiIlFUIIuIiIiIRFGBLCIiIiISRQWyiIiIiEgUFcgiIiKemFkXM/vRzOaZ2e0x1lc1s/Hh+n+bWcuodXeEy380s/OL69PMDgv7mBf2WWW/g48bBy1bQlpacDtu3H53JVIRqUAWEZGUtm3bNlasWFHu2zWzdOCfwAXA0cB1ZnZ0gWY3Ahucc0cCfweGh889GugGHAN0AUaYWXoxfQ4H/h72tSHsu/TGjYMBA2DRInAuuB0wQEWyJBUVyCKlNHLGSBZsWOA7hojEyYsvvkjz5s354YcfynvTHYB5zrkFzrlc4EXg0gJtLgWeC+9PAM41MwuXv+ic2+Wc+xmYF/YXs8/wOeeEfRD2edl+pb7rLti+/ZfLtm8PloskCRXIIqXww9ofyJqcxWs/vOY7iojESXZ2NkceeSRt2rQp7003A5ZEPV4aLovZxjm3B9gEHFTEcwtbfhCwMeyjsG0BYGYDzGymmc1cs2bNvg0WL479agpbLpKAVCCLlELu3lwubn0x1x17ne8oIhIHCxcuZNq0afTq1YvgIKs450Y75zKdc5mNGsW4Cm9GRuwnFrZcJAGpQBYphXaN2zHpukk0rd3UdxQRiYNx4bjZnj17+tj8MqBF1OPm4bKYbcysElAXWFfEcwtbvg6oF/ZR2LZKZtgwqFHjl8tq1AiWiyQJFcgiJfTzhp9ZvElfIYokC+cc2dnZnHHGGbRs2dJHhBlAq3B2iSoEJ91NKtBmEtAnvH8V8IFzzoXLu4WzXBwGtAK+LKzP8Dkfhn0Q9vn6fqXu0QNGj4ZDDwWz4Hb06GC5SJJQgSxSQsM+GcZxI49j155dvqOISBzMmDGDn376id69e3vZfjge+LfAFGAO8JJz7jszG2pml4TNngEOMrN5wGDg9vC53wEvAd8D7wA3Oef2FtZn2NdtwOCwr4PCvvdPjx6wcCHk5QW3Ko4lyVQqvomI7Ni9g5e/f5kr2l5B1UpVfccRkTiIRCJUq1aNq666qvjGZcQ5NxmYXGDZvVH3dwJXF/LcYcA+4xpi9RkuX0Awy4WIFENHkEVKYNKPk9i8azO92vXyHUVE4iA3N5cXXniBSy65hLp16/qOIyIVjApkkRLIzsmmeZ3mnNXyLN9RRCQO3nnnHdatW+dteIWIVGwqkEWKsXnXZj74+QN6HteTNNMuI5IMIpEIjRo1onPnzr6jiEgFpDHIIsWoU7UOC3+/UHOkiiSJDRs2MGnSJAYOHEjlypV9xxGRCkgFskgJNK7V2HcEEYmTl19+mdzcXA2vEJFC6ftikSJ8t/o7znnuHL5f873vKCISJ5FIhLZt23LiiSf6jiIiFZQKZJEiRHIiTFs0jYY1GvqOIiJxsGDBAj799FNdWlpEiqQCWaQQeS6Pcf8ZR5cju3BwzYN9xxGROBg7dixmRg9d2EJEiqACWaQQHy38iKWbl2ruY5Ek4ZwjEolw1llnkZGR4TuOiFRgKpBFChHJiVCnah0uaXNJ8Y1FpMKbPn068+bNo1cvfegVkaJpFguRQvy6+a9p3aA11StX9x1FROIgEolQvXp1rrzySt9RRKSCU4EsUogBvxrgO4KIxMmuXbsYP348l112GXXq1PEdR0QqOA2xEIlhyrwpbNm1xXcMEYmTyZMns379eg2vEJESUYEsUsCKLSu48PkLGf7ZcN9RRCROIpEIjRs35rzzzvMdRUQSgApkkQKe/8/z5Lk8erbr6TuKiMTB+vXrefPNN+nevTuVKmlkoYgUTwWySAGRnAgnHXISRzU8yncUEYmD8ePHs3v3bg2vEJESK7ZANrMxZrbazL4tZH0PM8sxs/+Y2edmdnz8Y4qUj5xVOXyz6ht6H9/bd5Qypf1aUkkkEuHYY4+lffv2vqOISIIoyRHkZ4EuRaz/GTjTOXcc8AAwOg65RLyYOn8qldIq0e3Ybr6jlLVn0X4tKWDevHl88cUXurS0iJRKsQWyc24asL6I9Z875zaED6cDzeOUTaTcDTllCPNvnk/DGg19RylT2q8lVUQiEcyM7t27+44iIgkk3mOQbwTeLmylmQ0ws5lmNnPNmjVx3rRIfGTU1SVoC9B+LQnJOcfYsWM599xzad5cn/FEpOTiViCb2dkEb6S3FdbGOTfaOZfpnMts1KhRvDYtEhcD3xzITW/d5DtGhaL9WhLZ559/zoIFC3RynoiUWlwKZDNrBzwNXOqcWxePPkXK09bcrURyIuzO2+07SoWh/VoSXXZ2NjVq1OCKK67wHUVEEswBF8hmlgFMBHo553468Egi5e/VOa+yffd2erXTkSbQfi2Jb+fOnbz00ktcccUV1KpVy3ccEUkwxc6YbmYvAGcBDc1sKXAfUBnAOTcKuBc4CBgRniG8xzmXWVaBRcpCdk42h9U7jFMzTvUdpVxov5Zk99Zbb7Fx40YNrxCR/VJsgeycu66Y9f2AfnFLJFLOlm1exvsL3ufuM+4mzVLj2jnaryXZZWdn07RpU84991zfUWIyswbAeKAlsBC4JmrmmOh2fYC7w4cPOueeC5f/imC6xurAZOD3zjlXWL9mdhTwL+BE4C7n3CNl9dpEkkFqVAMiRUizNG455ZakvziISKpYu3YtkydPpkePHqSnp/uOU5jbgfedc62A98PHvxAWu/cBHYEOwH1mVj9cPRLoD7QKf/LnNS+s3/XAzYAKY5ESUIEsKa9p7ab89by/cmSDI31HEZE4GD9+PHv27KnowysuBZ4L7z8HXBajzfnAVOfc+vDo8lSgi5k1Beo456Y75xyQHfX8mP0651Y752YAOhNZpARUIEtKm7d+Hu/Of5e9eXt9RxGROMnOzqZdu3a0a9fOd5SiNHbOrQjvrwQax2jTDFgS9XhpuKxZeL/g8pL2KyLFUIEsKW3EjBF0fb4rG3du9B1FROLgxx9/5Msvv6R3b/9DpszsPTP7NsbPpdHtwqPALt7b399+dfEfkRKcpCeSrPbk7eH5/zxP19ZdOajGQb7jiEgcjB07lrS0tApxaWnnXKfC1pnZKjNr6pxbEQ6ZWB2j2TKC2WbyNQc+Cpc3L7B8WXi/JP0Wl3s0MBogMzMz7oW7SCLQEWRJWe8teI9V21Zp7mORJJGXl0ckEuG8886jadOmvuMUZxLQJ7zfB3g9RpspQGczqx+enNcZmBIOodhsZidbMA9j76jnl6RfESmGCmRJWdnfZFO/Wn0ubHWh7ygiEgeffvopixYtqugn5+V7GDjPzOYCncLHmFmmmT0N4JxbDzwAzAh/hobLALIIrnQ5D5gPvF1Mv03COc8HA3eb2VIzq1P2L1MkMWmIhaSkvXl7mbF8Bt2O7UbVSlV9xxGROIhEItSsWZPLLos1IUTFEl6+fZ9Jmp1zM4mag9w5NwYYU0i7Y0vR70p+OSxDRIqgAllSUnpaOj/c9ANbc7f6jiIicbBjxw5eeuklrrrqKmrWrOk7jogkOA2xkJSU5/JIT0unbrW6vqOISBy88cYbbN68OVGGV4hIBacCWVLO4k2LOfQfhzJ1/lTfUUQkTiKRCM2aNeOss87yHUVEkoAKZEk543LGsXTzUo5ocITvKCISB6tXr+btt9+mZ8+eFfnS0iKSQFQgS0pxzhHJiXBqi1M5vP7hvuOISBy8+OKL7N27V8MrRCRuVCBLSvlqxVfMWTuH3sf7v8qWiMRHJBLhhBNO4JhjjvEdRUSShApkSSnZ32RTJb0KVx99te8oIhIHc+bMYebMmRXi0tIikjw0zZuklCvaXsERDY6gfvX6vqOISBxEIhHS09O57rrrfEcRkSSiAllSypktz+TMlmf6jiEicZCXl8fYsWPp3LkzjRs39h1HRJKIhlhIyhibM5bvVn/nO4aIxMnHH3/MkiVLNLxCROJOBbKkhE07N9FvUj9GzhzpO4qIxEkkEqF27dpceumlvqOISJJRgSwpYcL3E9i1dxe92mkaKJFksH37diZMmMBVV11F9erVfccRkSSjAllSQnZONq0Pak2HZh18RxGROHj99dfZsmWLhleISJlQgSxJb+HGhUxbNI1e7XphZr7jiEgcRCIRMjIyOOOMM3xHEZEkpAJZkt7slbOpWbkmPdv19B1FROJg5cqVTJkyhR49epCWprcxEYk/TfMmSe+yoy5jzZ/WUL2yximKJIMXXniBvLw8XVpaRMqMPnpLUsvdmwug4lgkiUQiETIzM2nbtq3vKCKSpFQgS1IbMmUIp445lTyX5zuKiMTBt99+y9dff62jxyJSpootkM1sjJmtNrNvC1l/lJl9YWa7zOyW+EcU2T+5e3N58bsXaV6nOWmmz4LRtF9LoopEIlSqVIlu3br5jiIiSawkVcOzQJci1q8HbgYeiUcgkXh5Z947rN2+lt7tNA1UDM+i/VoSzN69exk3bhxdunTh4IMP9h1HRJJYsQWyc24awZtlYetXO+dmALvjGUzkQEVyIjSq0YjOR3T2HaXC0X4tiejDDz9k2bJlGl4hImWuXL93NrMBZjbTzGauWbOmPDctKWbDjg1M+nES1x17HZXTK/uOk9S0X0t5iUQi1K1bl4svvth3FBFJcuVaIDvnRjvnMp1zmY0aNSrPTUuKqVapGk92fZLfZP7Gd5Skp/1aysO2bdt45ZVXuPrqq3VpaREpc5oHWZJS9crV6du+r+8YIhInr776Ktu2bdPwChEpFzq1X5LOwo0L+cf0f7Bx50bfUUQkTiKRCC1btuS0007zHUVEUkBJpnl7AfgCaGNmS83sRjMbaGYDw/VNzGwpMBi4O2xTp2xjixTuudnPMXjKYLbs2uI7SoWl/VoSyfLly3nvvffo2bNn0lxa2swamNlUM5sb3tYvpF2fsM1cM+sTtfxXZvYfM5tnZo+bmRXVr5n1MLOc8Dmfm9nx5fNKRRJTsUMsnHPXFbN+JdA8bolEDoBzjkhOhLNankWLui18x6mwtF9LInn++eeT8dLStwPvO+ceNrPbw8e3RTcwswbAfUAm4IBZZjbJObcBGAn0B/4NTCaYtvHtIvr9GTjTObfBzC4ARgMdy+F1iiSk5PgoLhKavnQ68zfMp/fxmvtYJFlEIhE6duxI69atfUeJp0uB58L7zwGXxWhzPjDVObc+LIqnAl3MrClQxzk33TnngOyo58fs1zn3edgHwHT0AVikSCqQJalEciJUr1SdK9te6TuKiMTBN998Q05OTrIdPQZo7JxbEd5fCTSO0aYZsCTq8dJwWbPwfsHlJe33RoKjzTFp6kYRzWIhSWb1ttVc0fYKalet7TuKiMRB/qWlr732Wt9RSs3M3gOaxFh1V/QD55wzMxfv7cfq18zOJiiQCz3b0Tk3mmAIBpmZmXHPJZIIVCBLUplwzQT25u31HUNE4mDv3r08//zzXHTRRTRs2NB3nFJzznUqbJ2ZrTKzps65FeGQidUxmi0Dzop63Bz4KFzevMDyZeH9Qvs1s3bA08AFzrl1+/GSRFKGhlhI0siftSI9Ld1zEhGJh/fff58VK1Yk4/AKgElA/qwUfYDXY7SZAnQ2s/rhbBSdgSnhEIrNZnZyOHtF76jnx+zXzDKAiUAv59xPZfGCRJKJCmRJCuu2r6PJo014atZTvqOISJxkZ2dTr149unbt6jtKWXgYOM/M5gKdwseYWaaZPQ3gnFsPPADMCH+GhssAsgiOBs8D5vO/McUx+wXuBcdvcaMAACAASURBVA4CRpjZbDObWcavTyShaYiFJIWXvnuJ7bu3c1Kzk3xHEZE42Lp1K6+++iq9evWiatWqvuPEXTjE4dwYy2cC/aIejwHGFNLu2FL02y+6XxEpmo4gS1LIzsnm2IOP5fjGmvteJBlMnDiR7du3J+vwChGp4FQgS8Kbu24u05dOp1e7XoQXkxKRBJednc3hhx/OKaec4juKiKQgFciS8MbmjMUwehzXw3cUEYmDpUuX8sEHH9Crlz70iogfGoMsCe/GE2+k1UGtaFanWfGNRaTCe/7553HO0bNnT99RRCRF6QiyJLyMuhn0bKc3UpFk4JwjOzubU045hSOPPNJ3HBFJUSqQJaGNmDGC1354zXcMEYmT2bNn89133+nkPBHxSgWyJKyde3Zy5/t38sqcV3xHEZE4iUQiVKlShWuuucZ3FBFJYSqQJWG9+dObbNq1iV7tdKRJJBns2bOH559/nq5du9KgQQPfcUQkhalAloQVyYnQtFZTzj1snznxRSQBTZ06lVWrVml4hYh4pwJZEtKabWuYPHcyPY7rQXpauu84IhIHkUiEBg0acOGFF/qOIiIpTgWyJKRFmxZxRP0j6HW8jjSJJIPNmzfz6quv0q1bN6pUqeI7joikOM2DLAkp85BM5tw0RxcREEkSr7zyCjt37tTwChGpEHQEWRLOpp2b2Llnp4pjkSQSiURo1aoVHTt29B1FREQFsiSev372V5o/1pxtudt8RxGROFi8eDEffvihLi0tIhWGCmRJKHkuj7H/GUvmIZnUrFLTdxwRiYNx48YB6NLSIlJhqECWhPLJok9YvGkxvY/v7TuKiMSBc45IJMJpp53GYYcd5juOiAigAlkSTCQnQq0qtbjsqMt8RxGROJg1axZz5syhd2996BWRikMFsiSMHbt38PL3L3PV0VdRo3IN33FEJA4ikQhVq1bl6quv9h1FROS/NM2bJIxqlaoxtddUalep7TuKiMTB7t27eeGFF7j44oupV6+e7zgiIv9V7BFkMxtjZqvN7NtC1puZPW5m88wsx8xOjH9METAzOjTrQNtGbX1HSXjar6UimDJlCmvWrNHwChGpcEoyxOJZoEsR6y8AWoU/A4CRBx5L5JdWb1vNgDcGMH/9fN9RksWzaL8WzyKRCA0bNqRLl6L+FEVEyl+xBbJzbhqwvogmlwLZLjAdqGdmTeMVUATghf+8wFNfPcXOPTt9R0kKPvbrWbNm8fHHHx9IF5JENm7cyOuvv063bt2oXLmy7zgiIr8QjzHIzYAlUY+XhstW7G+Hu/fu5tnZz3JovUPpfETnA80nSSA7J5sTm57IMQcf4ztKqijxfm1mAwiOMpORkRGzM+ccN9xwA2lpaXz11Ve6GIQwYcIEdu3apeEVIlIhlessFmY2wMxmmtnMNWvWFNouPS2d4Z8N58FpD5ZjOqmovlv9HV+t+Ipe7Xr5jiIxOOdGO+cynXOZjRo1itnGzMjKymL27NlMnz69nBNKRRSJRGjTpg2ZmZm+o4iI7CMeBfIyoEXU4+bhsn2U5I0UIM3SGJg5kE8Wf8K3q2OeQyQpJJITId3Sue7Y63xHSSUl3q9LqkePHtSuXZuRIzWcOdUtXLiQadOm6dLSIlJhxaNAngT0Ds96PxnY5Jzb7+EV+a5vfz1V06sycobeTFNdjco16H5cdxrXauw7SiqJ+35dq1Ytevfuzfjx41m7dm18UkpCGjt2LJDal5Y2swZmNtXM5oa39Qtp1ydsM9fM+kQt/5WZ/SecaeZxCz9pFNavmV0azkgzO/wm97TyeaUiiakk07y9AHwBtDGzpWZ2o5kNNLOBYZPJwAJgHvAUkBWPYAfVOIhrj72W7JxstuzaEo8uJUHde+a9ZF+e7TtGUvG1Xw8aNIjc3FzGjBkTj+4kAeVfWvrMM8/k0EMP9R3Hp9uB951zrYD3w8e/YGYNgPuAjkAH4L6oQnok0J//zTaTPxVIYf2+DxzvnGsP3AA8XRYvSiRZFHuSnnOuyO+1nXMOuCluiaJkZWYxZ80clm1ZxlFVjyqLTUgFN2/9PI6of4S+ho0zX/v1McccwxlnnMGoUaO45ZZbSEvTxTxTzYwZM/jpp5+49dZbfUfx7VLgrPD+c8BHwG0F2pwPTHXOrQcws6lAFzP7CKgTzjCDmWUDlwFvF9avc25rVL81ARfPFyOSbCr0u1PH5h35sv+XHNVQxXEq2r57Oyc8eQK3Tk35N9KkkpWVxc8//8yUKVN8RxEPsrOzqVatGldddZXvKL41jhq2tBKINYassNlkmoX3Cy4vsl8zu9zMfgDeIjiKHFNJT6gXSWYVukDOt2HHBpZsWlJ8Q0kqr/3wGltzt3JR64t8R5E4uvzyy2ncuLFO1ktBubm5vPjii1x66aXUrVvXd5wyZ2bvmdm3MX4ujW4XfmMT9yO6Bft1zr3qnDuK4GjzA0U8r0Qn1IskswpfIO/N28sxI47htvcKfvMkyS6SEyGjbgZnHHqG7ygSR1WqVKFfv368+eabLFq0yHccKUfvvPMO69atS5m5j51znZxzx8b4eR1YlX/xnfB2dYwuCptNZll4v+ByStJveKGgw82s4QG+RJGkVeEL5PS0dK455homfD+BVVtX+Y4j5WTFlhW8O/9deh7XkzSr8H+mUkoDBgzAzBg9erTvKFKOsrOzOfjgg+ncWReAIpgpJn9Wij7A6zHaTAE6m1n98OS8zsCUcAjFZjM7OZy9onfU82P2a2ZHRs10cSJQFVgX/5clkhwSovIYmDmQ3Xm7GfO1znxPFS999xJ5Lo9ex+viIMkoIyODrl278vTTT5Obm+s7jpSDDRs28MYbb3DddddRqVI8LuKa8B4GzjOzuUCn8DFmlmlmTwOEJ+c9AMwIf4bmn7BHMLPM0wQzzcwnOEGv0H6BK4FvzWw28E/g2nAIhojEkBD/Sx3V8CjOOewcRs0axa2n3kp6WrrvSFLGBp00iGMPPlYnaCaxQYMGMWnSJCZOnEi3bt18x5Ey9vLLL5Obm5sywyuK45xbB5wbY/lMoF/U4zHAPkeHwnbHlqLf4cDwA0stkjoS4ggyBFO+Ld60mC+WfuE7ipSDKulVOPfwff6PlyTSuXNnDj/8cEaMGOE7ipSD7Oxsjj76aE444QTfUUREipUwBfIlbS5hzk1zOC1DF/9Jdn//4u8MmzbMdwwpY2lpaQwcOJBPPvmEb7/VJeWT2fz58/nss890aWkRSRgJUyBXTq/836/bNWwqee3N28sjXzzCv5f923cUKQfXX389VatW1ZRvSW7s2LGYGT169PAdRUSkRBKmQAbIc3l0f6U7d39wt+8oUkY++PkDlm9ZTq92OjkvFTRs2JBrrrmGSCTCli26pHwyyr+09Nlnn02LFi2Kf4KISAWQUAVymqWRuzeXUbNGsWP3Dt9xpAxEciLUrVqXi9tc7DuKlJOsrCy2bNnCuHHjfEeRMjB9+nTmz59Pr1760CsiiSOhCmSArJOyWL9jPS9//7LvKBJnW3O3MnHORK455hqqVarmO46Uk44dO9K+fXtGjhyp4VNJKDs7m+rVq3PllVf6jiIiUmIJVyCf3fJs2hzUhpEzNWYx2WzYsYHOR3Smz/F9im8sScPMyMrKIicnh88//9x3HImjXbt2MX78eC6//HJq167tO46ISIklXIFsZgzKHMT0pdP5esXXvuNIHLWo24KJ107k1IxTfUeRcta9e3fq1Kmjk/WSTHZ2Nhs2bKBPH33oFZHEknAFMkCf9n24+/S7aVKrie8oEicf/vwh89bP8x1DPKlZsyZ9+vTh5ZdfZs2aNb7jSBxs2bKFe+65h1NPPZXzzjvPdxwRkVJJyAK5XrV6PHDOAzSt3dR3FImDbbnb6PVqL65//XrfUcSjQYMGkZuby5gxuqR8Mhg+fDirVq3iscce09zHIpJwErJAhmDqoEk/TmLSj5N8R5ED9OgXj7JsyzIePvdh31HEo7Zt23LWWWcxatQo9u7d6zuOHIAlS5bw6KOP0qNHDzp06OA7johIqSVsgWxmDPtkGLe9d5vOfE9gy7csZ/hnw7n66Ks19ljIyspi4cKFvPPOO76jyAG48847AXjooYc8JxER2T8JWyADZGVm8cPaH/ho4Ue+o8h+uueDe9iTt4eHO+noscBll11GkyZNdLJeApsxYwZjx45l8ODBZGRk+I4jIrJfErpAvuaYa2hQvYGmfEtQzjnqVavHLb++hcPrH+47jlQAlStXpn///kyePJmff/7ZdxwpJeccgwcP5uCDD+b222/3HUdEZL8ldIFcvXJ1bmh/A6/+8CrLtyz3HUdKycx49PxHGXbuMN9RpALp378/Zsbo0aN9R5FSmjhxIp9++ikPPvig5j0WkYSW0AUywMDMgWTUzeDnDTralEg+XvixhsZITC1atOCSSy7h6aefZteuXb7jSAnt2rWLW2+9leOOO44bbrjBdxwRkQOS8AXyEQ2OYN7v5ukErwSye+9uBrw5gEFvDWJvnmYrkH0NGjSItWvX8sorr/iOIiX0xBNPsGDBAh599FHS09N9xxEROSAJXyBD8FX9rj27WLJpie8oUgKjZo7ip3U/8ch5j5CepjdS2VenTp048sgjGTFihO8oUgJr167lgQce4MILL9RFQUQkKSRFgQxw6phTuWGSvtar6Dbs2MD9H99Pp8M7cWGrC33HkQoqLS2NgQMH8tlnn5GTk+M7jhRj6NChbN26lb/97W++o4iIxEXSFMiXH3U57y14j5/W/eQ7ihThwWkPsmHHBh457xFdXUuK1LdvX6pVq6Yp3yq4H374gREjRjBgwACOPvpo33FEROIiaQrkG0+8kcpplRk1c5TvKFKEw+sfzuBfD+b4Jsf7jiIV3EEHHcS1117L2LFj2bx5s+84Uohbb72VmjVrcv/99/uOIiISNyUqkM2si5n9aGbzzGyfyS3N7FAze9/McszsIzNrHv+oRWtSqwlXHn0l/5r9L7bv3l7em5cSuqnDTTzS+RHfMVJeIuzTEFxZb+vWrYwdO9bH5qUY77//Pm+88QZ33XUXBx98sO84IiJxU2yBbGbpwD+BC4CjgevMrOD3aI8A2c65dsBQ4C/xDloSgzIHsXHnRt748Q0fm5cifL7kc56d/Sx5Ls93lJSXSPv0SSedxIknnsjIkSN1SfkKZu/evQwZMoSWLVty8803+44jIhJXJTmC3AGY55xb4JzLBV4ELi3Q5mjgg/D+hzHWl4vTM05n+o3TueaYa3xsXgqR5/K4+e2buefDe9i5Z6fvOJJA+7SZkZWVxbfffsunn37qI4IU4rnnnuObb75h+PDhVKtWzXccEZG4KkmB3AyInj9tabgs2jfAFeH9y4HaZnZQwY7MbICZzTSzmWvWrNmfvEUyMzo276iTvyqYcTnjmLViFn859y/UqFzDdxyJ4z4NZb9fd+vWjbp16+pkvQpk69at3H333fz617/m6quv9h1HRCTu4nWS3i3AmWb2NXAmsAzY5woQzrnRzrlM51xmo0aN4rTpfd3+3u3c/La+8qsItu/ezp0f3EnmIZl0P6677zhSciXap6Hs9+uaNWvSt29fJkyYwKpVq+Lev5Te3/72N1asWMFjjz2mAxL7ycwamNlUM5sb3tYvpF2fsM1cM+sTtfxXZvaf8DyCxy38hyiuXzM7ycz2mNlVZfsKRRJbSQrkZUCLqMfNw2X/5Zxb7py7wjl3AnBXuGxj3FKW0qadm3jqq6dYt32drwgSeuyLx1i6eSmPdX6MNEuaSVMSXcLt04MGDWL37t2MGTPGVwQJLV26lL/97W9069aNk08+2XecRHY78L5zrhXwfvj4F8ysAXAf0JFgaNR9UQXvSKA/0Cr86VJcv+H5B8OBd8viBYkkk5JULDOAVmZ2mJlVAboBk6IbmFlDs/9WP3cAXt/FBp00iJ17dvLs7Gd9xhDghCYnMPjkwZx+6Om+o8j/JNw+3aZNG8455xxGjRrF3r26PLlPd911F3l5efzlL17O20wmlwLPhfefAy6L0eZ8YKpzbr1zbgMwFehiZk2BOs656S44ezU76vlF9fs74BVgdVxfiUgSKrZAds7tAX4LTAHmAC85574zs6FmdknY7CzgRzP7CWgMDCujvCXSrnE7Tss4jVGzRmnWBM8uan0Rj57/qO8YEiUR92kIpnxbvHgxb7/9tu8oKWvmzJlkZ2fzxz/+kZYtW/qOk+gaO+dWhPdXEuxnBRV2vkCz8H7B5YX2a2bNCM4nKHYwf1mfVyCSCEr0nbdzbrJzrrVz7gjn3LBw2b3OuUnh/QnOuVZhm37OuV1lGbokBmUOYt76eby34D3fUVLSNyu/4d4P72Vb7jbfUSSGRNynL7nkEpo2bcqIESN8R0lJzjmGDBlCo0aNuOOOO3zHSQhm9p6ZfRvj5xezwoRHgeM+j2GBfv8B3OZc8UeNyut8IZGKrJLvAGXlyrZX8u8O/+bQuof6jpJynHMMfncws1fO5o8n/5Ga1PQdSZJA5cqVGTBgAEOHDmXBggUcfvjhviOllNdee41p06YxcuRI6tSp4ztOQnDOdSpsnZmtMrOmzrkV4ZCJWMMelhF8m5OvOfBRuLx5geX55xEU1m8m8GJ4Ll9D4EIz2+Oce630r0wk+SXtWVNVK1Xl/y74P9o0bOM7Ssp5a+5bfPDzB9x/5v3Urx7zxGyR/dK/f3/S0tJ48sknfUdJKbm5udx6660cffTR9OvXz3ecZDEJyJ+Vog/weow2U4DOZlY/PDmvMzAlHEKx2cxODmev6B31/Jj9OucOc861dM61BCYAWSqORQqXtAVyvn8v/TcT50z0HSNl7N67m1vevYXWB7VmYOZA33EkyTRr1oxLL72UZ555hp07ddGZ8jJixAjmzZvHo48+SqVKSfvFY3l7GDjPzOYCncLHmFmmmT0N4JxbDzxAcGLtDGBouAwgC3gamAfMB94uql8RKZ2k/5/uwU8eZObymXRt3ZUq6VV8x0l6o2eN5sd1PzKp2yQqp1f2HUeS0KBBg5g4cSITJkygZ8+evuMkvfXr1zN06FDOP/98unTpUvwTpEScc+uAc2Msnwn0i3o8hhizyITtji1pvwXa9C19YpHUkvRHkAdlDmLl1pW89oO+SSoPp2acyq2n3ErX1l19R5Ekdc4559C6dWudrFdOhg4dyqZNm3jkkUd8RxERKTdJXyCff8T5HFbvMEbM0JtpeWjfpD3Dzxuuq2tJmUlLS2PgwIF88cUXzJ4923ecpPbTTz/xz3/+k/79+3PssfscrBQRqZjGjYOWLSEtLbgdN67UXSR9gZyels7AzIF8vOhjvl/zve84SWvBhgVc//r1rNy60ncUSQF9+/alevXqjBxZ7JSucgBuvfVWqlWrxp///GffUURESmbcOBgwABYtAueC2wEDSl0kJ32BDHDDCTfQtFZTflz7o+8oSeu2927jpe9e0oVZpFzUr1+fbt26MW7cODZt2uQ7TlL66KOPeP3117nzzjtp3DjWNSxERCqgu+6C7dt/uWz79mB5KaREgdywRkOW/HEJl7e93HeUpPTp4k+Z8P0Ebjv1Ng6pfYjvOJIisrKy2LZtG5FIxHeUpJOXl8fgwYPJyMjgD3/4g+84IiIlt3hx6ZYXIiUKZAiGWuS5PFZsWVF8YymxPJfHkHeH0Kx2M4b8eojvOJJCMjMzyczMZOTIkQQXDJN4iUQifP311zz88MNUr17ddxwRkZLLyCjd8kKkTIEMcPXLV9NlXBe9mcbRi9++yJfLvuShcx+iZhVdMU/KV1ZWFt9//z3Tpk3zHSVpbNu2jTvvvJOOHTvSrVs333FEREpn2DCoUeOXy2rUCJaXQkoVyOcfcT45q3L4YukXvqMkjbNbns39Z95Pz3aaj1bK37XXXkv9+vV1sl4cPfLIIyxfvpzHHntMs9GISOLp0QNGj4ZDDwWz4Hb06GB5KaRUgdz9uO7UqVpHU77FUdPaTbnvrPtIs5T6U5IKokaNGvTt25dXXnmFlSs1g8qBWrZsGX/961+5+uqrOeWUU3zHERHZPz16wMKFkJcX3JayOIYUK5BrValFn+P78PL3L7Nm2xrfcRLayq0ruWDcBZo6T7wbOHAge/bs4ZlnnvEdJeHdfffd7Nmzh4cf1tWJRSS1pVSBDDAwcyC5e3MZmzPWd5SEdu+H9/LegveonKbLSYtfrVu3plOnTjz55JPs3bvXd5yE9fXXX/Pcc8/x+9//nsMPP9x3HBERr1KuQD660dF80PsDftvht76jJKycVTk88/Uz/Pak39LqoFa+44iQlZXFkiVLeOutt3xHSUjOOYYMGcJBBx3EnXfe6TuOiIh3KVcgA5x92NlUTteRz/3hnGPIu0OoW7Uu95x5j+84IgBcfPHFNGvWjBEjdH7B/njjjTf48MMP+fOf/0y9evV8xxER8S4lC2SAJ758gv6T+vuOkXCmzJ/Cewve494z76VB9Qa+44gAUKlSJQYMGMCUKVOYP3++7zgJJTc3l1tuuYW2bdsyYMAA33FERCqElC2Q12xbwzNfP8PPG372HSWhnHHoGfzj/H+QdVKW7ygiv9CvXz/S09MZNWqU7ygJZdSoUcydO5dHHnmESpUq+Y4jIlIhpGyB3P9X/UmzNJ6c9aTvKAmlRuUa/P7k31MlvYrvKCK/cMghh3D55ZczZswYduzY4TtOQli/fj33338/nTp14oILLvAdR0SkwkjZArl5neZc0uYSnvn6GXbt2eU7ToW3aecmTn76ZD5a+JHvKCKFGjRoEOvXr+fll1/2HSUhPPjgg2zcuJFHH31UFwUREYmSsgUyQNZJWazdvpYJ30/wHaXCe+iTh/hy2ZfUqVrHdxSRQp199tm0adNGV9YrgXnz5vHEE09w44030q5dO99xREQqlJQukM857BwG/mogRzQ4wneUCu3nDT/zj3//g97H9+bEpif6jiNSKDNj0KBBTJ8+na+++sp3nArttttuo0qVKjzwwAO+o4iUnXHjoGVLSEsLbseN851IEkRKF8hplsbIriM5ufnJvqNUaHe8fwfpls6wc4b5jiJSrD59+lC9enUdRS7CtGnTmDhxInfccQdNmjTxHUekbIwbBwMGwKJF4FxwO2CAimQpkZQukPPNWz9PwywK8dWKrxj/3Xj+dMqfaFanme84IsWqV68e3bt35/nnn2fjxo2+41Q4eXl5DB48mBYtWjB48GDfcUTKzl13wfbtv1y2fXuwXKQYKpAJxtf2fa0vm3Zu8h2lwmnfpD0vXPkCfzr1T76jiJRYVlYW27dvJzs723eUCmfcuHHMmjWLhx56iOrVq/uOI1J2Fi8u3XKRKCUqkM2si5n9aGbzzOz2GOszzOxDM/vazHLM7ML4Ry07WSdlsW33NiI5Ed9RKpQ8l0eapdHt2G7UqlLLdxyJo2Tfp0888UQ6dOjAyJEjcc75jlNhbN++nTvuuIPMzEy6d+/uO45I2crIKN1ykSjFFshmlg78E7gAOBq4zsyOLtDsbuAl59wJQDcgoa73mnlIJicdchIjZ+rNNN+O3Ts44ckTGJsz1ncUibNU2KchOIr8ww8/8NFHH/mOUmE89thjLFu2jMcee4y0NH2BKElu2DCoUeOXy2rUCJaLFKMk/0N2AOY55xY453KBF4FLC7RxQP78X3WB5fGLWD6yTsri+zXfM23RNN9RKoT/+/f/kbMqh0NqH+I7isRfSuzT11xzDQ0aNNDJeqEVK1bw8MMPc+WVV3L66af7jiNS9nr0gNGj4dBDwSy4HT06WC5SjJIUyM2AJVGPl4bLot0P9DSzpcBk4HexOjKzAWY208xmrlmzZj/ilp1rj7mWJrWa8O3qb31H8W71ttU89MlDXNLmEs457BzfcST+4rZPQ8Xdr6tXr87111/Pq6++yooVK3zH8e6ee+4hNzeX4cOH+44igJk1MLOpZjY3vK1fSLs+YZu5ZtYnavmvzOw/4TCpxy280kth/ZrZWWa2ycxmhz/3ls8r9axHD1i4EPLyglsVx1JC8fqO7TrgWedcc+BCIGJm+/TtnBvtnMt0zmU2atQoTpuOj+qVq7Pw9wu5qcNNvqN4d9+H97Fjzw7+2umvvqOIPyXap6Fi79cDBw5kz549PP30076jePXNN98wZswYbr75Zo44QvO+VxC3A+8751oB74ePf8HMGgD3AR0Jvvm5L6qQHgn0B1qFP11K0O8nzrn24c/QMnhNIkmjJAXyMqBF1OPm4bJoNwIvATjnvgCqAQ3jEbA8Va1UFYD1O9Z7TuLP0s1LeeqrpxiUOYg2Ddv4jiNlI2X26SOPPJLOnTvz5JNPsmfPHt9xvHDOMXjwYOrXr89dmt6qIrkUeC68/xxwWYw25wNTnXPrnXMbgKlAFzNrCtRxzk13wYkz2VHPL0m/IlKMkhTIM4BWZnaYmVUhOGFnUoE2i4FzAcysLcGbacX5rrUUhkwZwvGjjmdPXmq+mTav05yP+37MvWemxrdvKSql9umsrCyWLVvGm2++6TuKF2+99RYffPAB999/P/Xrx/wWX/xo7JzLH/uzEmgco01hw6GahfcLLi+u31+b2Tdm9raZHVNYsIo6bEqkPBVbIDvn9gC/BaYAcwjObP/OzIaa2SVhsyFAfzP7BngB6OsSdDqI0w89naWbl/LWT2/5jlLucvfmAnBqxqk0rJFwBwulhFJtn77oooto3rw5I0Yk3EQcB2z37t3ccssttG7dmoEDB/qOk3LM7D0z+zbGzy9Oig33rbjvXwX6/Qo41Dl3PPD/gNeKeF6FHTYlUl4qlaSRc24ywYk60cvujbr/PXBqfKP50bV1V5rXac6ImSO49KiCJ/Ynrz15e+jwVAeuPvpq7jpDX8Mmu1TapytVqsRvfvMb7rnnHubOnUurVq18Ryo3o0eP5scff2TSpElUrlzZd5yU45zrVNg6M1tlZk2dcyvCjqD9fwAAEnBJREFUIROrYzRbBpwV9bg58FG4vHmB5fnDpGL265zbHJVrspmNMLOGzrm1+/HSRJKeJsIsoFJaJQacOIB357/L3HVzfccpN8989QzfrPqGto3a+o4iEnf9+vWjUqVKjBo1yneUcrNx40buu+8+zjnnHLp27eo7juxrEpA/K0Uf4PUYbaYAnc2sfnhyXmdgSjiEYrOZnRzOXtE76vkx+zWzJlEzXXQgeP9fF/+XJZIcVCDH0O/EflRKq8RTXz3lO0q52LxrM/d8eA+nZ5zO5Udd7juOSNw1adKEK664gn/961/s2LHDd5xyMWzYMNavX8+jjz5KWBdJxfIwcJ6ZzQU6hY8xs0wzexrAObceeIDgvIEZwNBwGUAW8DQwD5gPvF1Uv8BVwLfhsKnHgW6JOmxKpDyUaIhFqmlauylvXvcmp7Q4xXeUcvGXT/7Cmu1rmHz+ZL2RStIaNGgQL730EuPHj6dv376+45Sp+fPn8/jjj9O3b1/at2/vO47E4JxbR3gibIHlM4F+UY/HAGMKaXdsKfp9AnjiwFKLpA4dQS7E+UeeT+2qtX3HKHNbc7cycuZIerXrReYhmb7jiJSZM888k7Zt26bElfVuv/12KlWqxIMPPug7iohIQlKBXIRXvn+FPq/1Kb5hAqtVpRZf/+ZrhnfS1bUkuZkZgwYN4ssvv2TWrFm+45SZTz/9lAkTJnDbbbdxyCG6VLyIxNG4cdCyJaSlBbfjxvlOVGZUIBdhxdYVZH+TzczlM31HKRNbdm0B4LD6h9G0dlPPaUTKXu/evalRo0bSHkXOy8tjyJAhNPv/7d17cFR1msbx70u4OAk3CXK/CBIVEUhiYHXUDY6WImsJljKFG7wASgyyOlvsjtRQCqPFuIUOuM4KGtEFd2Eyo2sprFdm3BmnXEEwAUFZUKYQwSAIeOMihLz7Rx9iExLohE6f7s7zqepK9+lfd54+6Tfnzcmvz+nZk2nTpoUdR0TSyZIlMHkyfPopuEe+Tp6ctk2yGuSTuGXILWS2ymTB6vTbmLo71/znNUx4eULYUUQSpkOHDhQVFbF06VL27dsXdpy4Kysr47333uNXv/oVWVlZYccRkXQyYwYcOHD8sgMHIsvTkBrkk+hwRgfGDx7P0g1L2XcwvTamz3/0PO9uf5fLel8WdhSRhJoyZQoHDx5k8eLFpx6cQg4ePMj06dPJz89n/PjxYccRkXSzbVvDlqc4NcinUDKshENVh1i0dlHYUeLmUNUh7vvDfQzpOoTbc28PO45IQuXm5nLJJZewYMEC0ukoV/PmzeOzzz5j7ty5tGihX+0iEmd9+jRseYrTb9FTyO2WS0lBCTnZ6XP2rcdXPc7Wr7by66t/TUaLjLDjiCRcSUkJmzdv5q233go7Slzs3LmThx9+mDFjxlBYWBh2HBFJR7NnQ2bm8csyMyPL05Aa5BjM/7v5XHduepyJ6mj1URaWL+S6c6/jqv71ngVVJK2NHTuW7OzstPmw3gMPPMChQ4eYM2dO2FFEJF0VFUFpKfTtC2aRr6WlkeVpSCcKidGXB77kL5/+hRsGpvaZ5jJaZLBm8pqaI1iINEdnnHEGEydOZO7cuezYsYOePXuGHanR1q9fzzPPPMM999xDTk76/KdLRJJQUVHaNsS1aQ9yjB7930cZ+/xYtn+zPewojfblgS+pqq6ifZv29Gyfug2BSDwUFxdTXV3NwoULw47SaO7OtGnT6NChA/fff3/YcURE0oYa5BhNvmgy1V7N0+8/HXaURrvtpdsYsWhEWn0wSaSxzjnnHK655hpKS0s5cuRI2HEa5fXXX2fFihXMnDmTTp06hR1HRCRtqEGOUf8z+3NtzrWUlpdy5GjqbUzf3PImr378KjecfwNmFnYckaQwZcoUPv/8c5YvXx52lAarqqpi2rRpDBgwgJKSkrDjiIikFTXIDTClYAo7v9vJS//3UthRGuRo9VGmvTmN/mf2Z+rwqWHHEUkao0aNok+fPsyfPz/sKA329NNPs3HjRh555BFat24ddhwRkbSiD+k1wMgBI+nXsR/vV77P2EFjw45Tr6rqKr469BWdMztz5OgRChcVsmHXBp4f+zxtWrYJO55I0sjIyKC4uJgZM2awadMmzjvvvLAj1Wvfvn2sXbuWiooKysvLWb58OYWFhYwePTrsaCIiaUcNcgNktMjgg5IPaNu6bdhRjlNeWc7qHaup2FlBeWU563et56r+V7H85uW0ymhFj3Y9eHDEg9w48Mawo4oknUmTJjFr1iyefPJJ5s2bF3Yc3J3KykrKy8upqKiouWzdurVmTI8ePRgxYgSPPvqopkyJiDQBNcgNdKw5/u7wdwlvlL/5/hvW7lxLeWU5+w7u45dX/BKAe167h3c+e4eOZ3Qkr1sedw+7m8v7XF7zuBd++kJCc4qkkq5du3LjjTeyaNEiZs+eTWbtA+E3oerqarZs2XJcI1xRUcGuXbtqxuTk5DB8+HCKi4vJy8sjLy+PLl26JCyjiEhzpAa5ER5b+RgPvf0Q2362jazWWU3yPXbv381ZWWcBMPfducxfPZ8t+7bU3H92x7OZOWImLawFT4x6gg5ndKBvh77amyTSCCUlJZSVlVFWVsbEiROb5HscOXKEjz76qKYJLi8vZ926dXz7beSY5C1btmTQoEGMGjWqphEeOnQo7du3b5I8IiJSPzXIjTCsxzD2HtzL0vVLufOiO0/7+Xbt38U7296pmSJRsbOCz7/9nC/+6Qu6ZHUhq1UWud1ymZA7gfzu+eR1z6Nb2241jx/abehpZxBpzi6//HIGDRrEggUL4tIg79+/n3Xr1h23V3jDhg0cPnwYgMzMTHJzc7n11ltrmuFBgwbRpo0+IyAikgzUIDfCj3v/mMFdBjN/zXzuyL8j5r22R6uPsnnP5ppG+K6CuxjQaQCvbH6Ficsm0sJaMLDzQH7S7yfkd8unZYvIj6e4oJjiguKmfEkizZqZUVJSwtSpU1m9ejXDhg2L+bF79uw5YYrEpk2bao43np2dTV5eHvfee29NM5yTk0NGRkZTvRwRETlNapAbwcyYMmwKJa+UsGrHKi7udfEJY76v+p7DRw/Trk07Nu7eyKRlk1j3xToOHDkAQJuMNhT2LWRApwFcd+51rJy0ksFdB5PZKnHzH0XkB7fccgv33XcfCxYsqLNBdne2b99+QjO8bdu2mjG9e/cmPz+fcePG1TTDvXr10tQnETm5JUtgxgzYtg369IHZs5vNKZ2TlRrkRioaXMTPV/yc+avnU9CjgJXbV1JRWUH5znIqKiv4cPeHPHTFQ0y/bDqdftSJVhmtuDP/zsgUiW55nN/5fFpltALgrKyzauYbi0g42rdvz/jx41m8eDFz5syp2TMcfTSJPXv2AJE/ks877zwuvfRSpk6dWtMMZ2dnh/wqRCTlLFkCkyfDgcgOND79NHIb1CSHSA1yI7Vr046ym8oY2nUo1V7Nlc9dyeGjh+mS1YW8bnmMyhlFYd9CALq27cqfb/9zyIlF5FRKSkp46qmn6N69O1VVVQC0bt2aCy+8kDFjxpCXl0d+fj5DhgwhK6tpPqArIs3MjBk/NMfHHDgQWa4GOTRqkE/DqJxRNdffHP8mOdk5dG/bXf9OFUlRQ4cOZdasWezdu7dmr/DAgQN1pjoRaTpR07RiWi4JoQY5TgrPLgw7gkjMzGwk8K9ABrDQ3f+l1v3zgCuCm5lAF3fvmNiU4Zg5c2bYEUSkOenTJzKtoq7lEpoWsQwys5FmtsnMPjGz6XXcP8/M1gaXzWb2Vfyjikg8mFkG8ARwLXABcLOZXRA9xt3/0d1z3T0X+A3wYuKTiog0A7NnQ+0TFGVmRpZLaE7ZIGtjKpJ2hgOfuPtf3f0wUAaMPsn4m4HfJiSZSDNhZp3MbIWZfRx8PbOecbcFYz42s9uill9kZuuDHVePWzC372TPa2Yjgh1ZH5qZPhiTLIqKoLQU+vYFs8jX0lLNPw5ZLHuQtTEVSS89gc+ibm8Plp3AzPoC/YC3EpBLpDmZDvzR3XOAPwa3j2NmnYCZwN8Q2RbPjGp4FwB3AjnBZeTJntfMOgLzgevdfRAwtolelzRGURFs3QrV1ZGvao5DF0uDHLeNqZlNNrM1ZrZm9+7dDc0qIok3DnjB3Y/WN0B1LdIoo4HFwfXFwJg6xlwDrHD3ve6+D1gBjDSz7kB7d1/pkTPSPBf1+Pqe9++BF919G4C774r3CxJJJzHNQW6Ak25M3b3U3QvcveCss3TcX5GQ7AB6R93uFSyryzhO8R8h1bVIo3R198rg+k6gax1j6ttB1TO4Xnv5yZ73XOBMM/uTmb1vZrfWF0x/9IrEdhSLhm5M7z7dUCLSpFYDOWbWj0gtjyOyd+k4ZnY+cCbwbmLjiaQHM/sD0K2Ou2ZE33B3NzOP9/ev9bwtgYuAK4EfAe+a2Up331zH40qBUoCCgoK45xJJBbE0yNqYiqQRd68ys6nAG0QO8/asu39oZg8Ca9x9WTB0HFAW/AtXRBrI3a+q7z4z+8LMurt7ZTBloq4pDzuAEVG3ewF/Cpb3qrX82I6r+p53O7DH3fcD+83sbWAocEKDLCIxTLFw9yrg2MZ0I/D7YxtTM7s+aqg2piIpwt1fdfdz3f0cd58dLHsgqjnG3We5+wkfHBKRuFgGHDsqxW3Ay3WMeQO42szODD6cdzXwRjCF4hszuzg4esWtUY+v73lfBi4zs5Zmlknkg38b4/2iRNKFhdXPmtluoI4jYx+nM/BlAuKcjlTICMoZb6fK2dfdm92E3BjqOl1+vslCOeMnloxxq2szywZ+D/QhUjM/dfe9ZlYA3OXudwTjJgK/CB42293/PVheACwiMl3iNeAfgikVdT5v8Jh/BiYA1UROEPRYDDlV04mlnPHV6G11aA1yLMxsjbsXhJ3jZFIhIyhnvKVKzmSTKutNOeMrFXKmQsZklCrrTTnjqznkjPdRLEREREREUpoaZBERERGRKMneIJeGHSAGqZARlDPeUiVnskmV9aac8ZUKOVMhYzJKlfWmnPGV9jmTeg6yiIiIiEiiJfseZBERERGRhFKDLCIiIiISJfQG2cxGmtkmM/vEzE44KYGZtTGz3wX3rzKzsxOfMqact5vZbjNbG1zuCCHjs2a2y8w21HO/mdnjwWv4wMzyE50xyHGqnCPM7OuodflAojMGOXqb2f+Y2Udm9qGZ3VvHmKRYp8lGdR3XjKrr+GVUTTeSajquGVXTcdRkde3uoV2InOZ2C9AfaA2sAy6oNWYK8GRwfRzwuyTNeTvwbyGvz78F8oEN9dw/isgB5Q24GFiVpDlHAP8d5roMcnQH8oPr7YickrX2zz0p1mkyXVTXcc+puo5fRtV049abajq+OVXT8c3ZJHUd9h7k4cAn7v5Xdz8MlAGja40ZDSwOrr8AXGlmlsCMEFvO0Ln728DekwwZDTznESuBjmbWPTHpfhBDzqTg7pXuXh5c/5bIaVl71hqWFOs0yaiu40h1HT+q6UZTTceRajq+mqquw26QewKfRd3ezokvqmaMu1cBXwPZCUlXR4ZAXTkBbgx23b9gZr0TE61BYn0dyeASM1tnZq+Z2aCwwwT/LswDVtW6K5XWaaKorhMrld6DSVPXqukGUU0nViq9B5OmpiG+dR12g5xOlgNnu/sQYAU//CUtDVdO5PzoQ4HfAC+FGcbM2gL/BfzM3b8JM4sknOo6fpKmrlXTzZpqOn6SpqYh/nUddoO8A4j+661XsKzOMWbWEugA7ElIujoyBE7I6e573P374OZC4KIEZWuIWNZ36Nz9G3f/Lrj+KtDKzDqHkcXMWhEpuCXu/mIdQ1JinSaY6jqxUuI9mCx1rZpuFNV0YqXEezBZahqapq7DbpBXAzlm1s/MWhOZ2L+s1phlwG3B9ZuAtzyYcZ1Ap8xZay7L9UTmwCSbZcCtwac5Lwa+dvfKsEPVZmbdjs1dM7PhRN6nif5FS5DhGWCju8+tZ1hKrNMEU10nVkq8B5OhrlXTjaaaTqyUeA8mQ00H37tJ6rplnHM2iLtXmdlU4A0inz591t0/NLMHgTXuvozIi/4PM/uEyGTxcUma8x4zux6oCnLenuicZvZbIp8q7Wxm24GZQKvgNTwJvErkk5yfAAeACYnOGGPOm4ASM6sCDgLjQvhFC3ApcAuw3szWBst+AfSJypoU6zSZqK7jS3UdV6rpRlBNx5dqOu6apK51qmkRERERkShhT7EQEREREUkqapBFRERERKKoQRYRERERiaIGWUREREQkihpkEREREZEoapBFRERERKKoQRYRERERifL/FOov+eeFHj4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "code", "metadata": { "id": "HGJqtjyijALS", "colab": { "base_uri": "https://localhost:8080/", "height": 203 }, "outputId": "2c767071-d58e-42af-e31b-834adaac6819" }, "source": [ "import pandas as pd\n", "d = {'time t_i': t, '4th Order Runge Kutta, w_i': w,'Exact':y,'Error |w-y|':np.round(np.abs(y-w),5)}\n", "df = pd.DataFrame(data=d)\n", "df" ], "execution_count": 7, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
time t_i4th Order Runge Kutta, w_iExactError |w-y|
00.01.0000001.0000000.00000
10.50.7135420.7130610.00048
21.00.7363420.7357590.00058
31.50.9467910.9462600.00053
42.01.2711001.2706710.00043
\n", "
" ], "text/plain": [ " time t_i 4th Order Runge Kutta, w_i Exact Error |w-y|\n", "0 0.0 1.000000 1.000000 0.00000\n", "1 0.5 0.713542 0.713061 0.00048\n", "2 1.0 0.736342 0.735759 0.00058\n", "3 1.5 0.946791 0.946260 0.00053\n", "4 2.0 1.271100 1.270671 0.00043" ] }, "metadata": {}, "execution_count": 7 } ] }, { "cell_type": "code", "metadata": { "id": "uK_MhzMfjG0V" }, "source": [ "" ], "execution_count": 7, "outputs": [] } ] }